Why the Z-Transform?

• The Fourier Transform assumes signals are stable and well-behaved
• But some biosignals or systems may not be absolutely summable
• The Z-Transform generalizes the Fourier Transform
• Useful for analyzing discrete-time systems, especially when stability and causality matter

Definition

Let \(x[n]\) be a discrete-time signal.

The Z-Transform is defined as:

\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]

Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)

Region of Convergence (ROC)

• The Z-Transform converges only for certain values of \(z\)
• The set of \(z\) for which the series converges is the ROC
• ROC is critical for system stability and causality

Z-Plane Representation

• Poles: values of \(z\) where \(X(z) \to \infty\)
• Zeros: values where \(X(z) = 0\)
• Visualization of poles and zeros helps in understanding system behavior

(-1.5, 1.5)

(-1.5, 1.5)

Relationship with Fourier Transform

If the ROC includes the unit circle, \(|z| = 1\), then:

\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]

So the Fourier Transform is a special case of the Z-Transform.

Example Signal

Let:

\[x[n] = a^n u[n]\]

Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))

Step 1: Z-Transform

Apply the definition:

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]

This is a geometric series:

\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]

What is a Geometric Series?

A geometric series is a sum of terms where each term is multiplied by the same constant:

\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]

The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.

Goal

Understand why this series:

\[\sum_{n=0}^{\infty} r^n\]

converges if and only if:

\[\boxed{|r| < 1}\]

Partial Sums

Let’s consider the sum up to the \(N\)-th term:

\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]

This has a closed-form expression:

\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]

Taking the Limit

To find the sum of the infinite series, take the limit as \(N \to \infty\):

\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]

Case 1: \(|r| < 1\)

If \(|r| < 1\), then:

\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]

So the sum becomes:

\[S = \frac{1}{1 - r}\]

✅ The geometric series converges.

Case 2: \(|r| \geq 1\)

• If \(r = 1\), then \(S_N = N + 1 \to \infty\)
• If \(r = -1\), the sum oscillates: \(1 - 1 + 1 - 1 + \cdots\)
• If \(|r| > 1\), then \(r^{N+1} \to \infty\)

❌ In all cases: the series diverges

Intuition

• When \(|r| < 1\), each term \(r^n\) gets smaller and smaller
• Their total sum settles to a finite number
• When \(|r| \geq 1\), the terms don’t vanish — the sum keeps growing or oscillating

Example: \(r = 0.5\)

\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]

Every term adds less. The sum “flattens out.”

Why This Matters

The Z-Transform often gives us geometric series like:

\[\sum_{n=0}^{\infty} (az^{-1})^n\]

This converges only if:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

So, understanding convergence of geometric series = understanding ROC in Z-transforms.

Summary

\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]

Step 2: Region of Convergence (ROC)

For convergence of the geometric series:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

Therefore, the ROC is:

\[\boxed{|z| > |a|}\]

Interpretation

• Causal signal (defined for \(n \geq 0\))
• ROC is outside the outermost pole
• Stable system only if ROC includes the unit circle \(|z| = 1\)

Example: \(a = 0.5\)

\[x[n] = (0.5)^n u[n]\]

Z-Transform:

\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]

Region of Convergence:

\[\boxed{|z| > 0.5}\]

Summary

Properties of the Z-Transform

• Linearity: \(a x[n] + b y[n] \to aX(z) + bY(z)\)
• Time shifting: \(x[n - k] \to z^{-k} X(z)\)
• Scaling in the z-domain: \(a^n x[n] \to X(z/a)\)
• Convolution: \(x[n] * h[n] \to X(z)H(z)\)

Example

Let \(x[n] = a^n u[n]\), where \(|a| < 1\)

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]

Difference Equations in DSP

A difference equation relates input and output values at different time steps.

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems

Z-Transform of Time-Shifted Terms

The Z-Transform turns time shifts into powers of \(z^{-1}\):

Step 1: Apply Z-Transform

Given:

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Apply \(\mathcal{Z} \{ \cdot \}\):

\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]

Step 2: Factor and Solve for \(H(z)\)

Group:

\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]

Divide both sides:

\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]

Example

Given:

\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]

Z-Transform:

\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]

Transfer Function:

\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]

Poles and Zeros

Let’s analyze \(H(z)\):

• Zeros: Roots of the numerator \(\Rightarrow z = 0.5\)
• Poles: Roots of the denominator \(\Rightarrow z = 0.9\)

Practice

Convert this equation:

\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]

Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane

Application in Biosignal Processing

• Analysis of digital filters for ECG, EEG, etc.
• Design of stable and causal filtering systems
• Useful in difference equation modeling of biosignals

Summary

• Z-Transform is a powerful tool for analyzing discrete systems
• Provides insight into stability, causality, and system behavior
• A generalization of the Fourier Transform
• Crucial in digital signal processing of biosignals
• Z-Transform converts difference equations into algebraic expressions
• Transfer function \(H(z)\) tells us how the system responds to inputs
• Key for digital filter design in biosignal processing

Next Steps

• Practice Z-Transform computations
• Pole-zero plotting exercises
• Application to real biosignal filtering problems

Filter Design

Filter Design

def design_filter(zeros=None, poles=None, gain=1.0):
    """
    Diseña un filtro digital a partir de ceros y/o polos y una ganancia.

    Parámetros:
    - zeros: lista de ceros (raíces del numerador), o None para no incluir
    - poles: lista de polos (raíces del denominador), o None para no incluir
    - gain: ganancia escalar del filtro

    Devuelve:
    - b: coeficientes del numerador
    - a: coeficientes del denominador
    """
    # Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
    if zeros:
        b = gain * np.poly(zeros)
    else:
        b = np.array([gain], dtype=float)

    # Si no se pasan polos, asumimos sistema FIR (a = [1])
    if poles:
        a = np.poly(poles)
    else:
        a = np.array([1.0], dtype=float)

    return b, a

Filter Design

def gaussian(x, mu, sigma, A):
    """
    Genera una función gaussiana.

    Parámetros:
    - x: array de tiempos
    - mu: posición central de la gaussiana
    - sigma: desviación estándar
    - A: amplitud
    """
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

Filter Design

def simulate_ecg(duration=10.0, fs=500, heart_rate=60):
    """
    Simula un ECG sintético basado en la superposición de ondas gaussianas.

    Parámetros:
    - duration: duración de la señal en segundos
    - fs: frecuencia de muestreo en Hz
    - heart_rate: frecuencia cardiaca en latidos por minuto

    Devuelve:
    - t: vector de tiempos
    - ecg: señal simulada de ECG en mV
    """
    dt = 1 / fs
    t = np.arange(0, duration, dt)
    rr = 60 / heart_rate  # intervalo RR en segundos

    # Inicializar señal
    ecg = np.zeros_like(t)

    # Parámetros de las ondas (posiciones relativas en segundos)
    # P wave
    p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
    # Q wave
    q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
    # R wave
    r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
    # S wave
    s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
    # T wave
    t_amp, t_dur, t_delay = 0.35, 0.142, 0.36

    # Generar cada latido
    for beat_start in np.arange(0, duration, rr):
        mask = (t >= beat_start) & (t < beat_start + rr)
        tb = t[mask] - beat_start
        ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
        ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
        ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
        ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
        ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)

    return t, ecg+0.1*np.cos(2*np.pi*60*t)

Filter Design

# Parámetros de simulación
DURATION = 10    # segundos
FS = 500         # Hz
HR = 70          # latidos por minuto

# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)

# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()

Filter Design

n = len(ecg_signal)
yf = np.fft.fft(ecg_signal)
xf = np.fft.fftfreq(n, 1/fs)[:n//2]
plt.plot(xf, 2.0/n * np.abs(yf[0:n//2]))
plt.grid()
plt.show()

Filter Design

fc = 60
b,a =design_filter(zeros=[np.exp(1j*2*np.pi*fc/fs), np.exp(-1j*2*np.pi*fc/fs)])
print(a)
print(b)

w, h = sig.freqz(a, b, worN=8000)
plt.plot(w/np.pi*fs/2, 20*np.log10(abs(h)))
plt.grid()
plt.xlabel('Frequency (Hz)')

Filter Design

Filter Design

You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),

\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]

By definition of the inverse DTFT,

\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]

Relating to sampling frequency

If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then

\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]

so you can write

\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]

where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).

Key takeaway

• \(h_d[n]\) is exactly the inverse‐DTFT of the ideal (“brick‐wall”) frequency specification.
• It produces a sinc-shaped impulse response of infinite length.
• Truncation (with a window) makes it realizable as a finite-length FIR.

Filter Design

The typical characteristics of classic windows are:

• Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.

• Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.

• Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).

• Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.

• Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).

• Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,

  \[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]

  and the normalized transition width Δω relates to the order M and β by

  \[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]

These formulas stem from Kaiser’s approximations and help in sizing the filter.

Windows Forms

(-100.0, 5.0)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
    boxcar,       # Rectangular
    bartlett,    # Triangular
    hann,        # Hann
    hamming,     # Hamming
    blackman,    # Blackman
    kaiser       # Kaiser
)

# Parameters
M = 64               # window length
beta = 8.6           # Kaiser parameter for moderate sidelobe attenuation
nfft = 512           # FFT length for frequency response
fs = 1.0             # normalized sampling rate

# Generate windows
windows = {
    'Rectangular': boxcar(M),
    'Bartlett':    bartlett(M),
    'Hann':        hann(M),
    'Hamming':     hamming(M),
    'Blackman':    blackman(M),
    f'Kaiser (β={beta})': kaiser(M, beta)
}

# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    # Compute normalized frequency response
    W = np.fft.fft(w, nfft)
    W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
    freqs = np.fft.fftfreq(nfft, d=1/fs)
    # Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
    idx = np.logical_and(freqs >= 0, freqs <= 0.5)
    plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

plt.show()

Filter Design

N = 101
fc = 30

# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)

# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w

# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)

# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)

# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()